En el presente documento se pretende ir comentando los pasos que se están siguiendo para conseguir predecir la cantidad de agua que llegará a la presa de belesar. Para ello se cuenta con la predicción meteorológica que diariamente se realiza con WRF.
Para la elaboracion del modelo predictivo se cuenta con un Histórico de datos 15 minutales aportados por la presa de Belesar en la que se aportan datos como Aportación, Nivel del embalse, Caudal turbinado, lluvia y algunas otras variables.
DHI<- readRDS(here::here('Data/Parques/Belesar/Historico/Historico_DHI_Belesar.RDS'))
DHI
Uno de los problemas que nos encontramos a la hora de tratar los datos es la presencia de valores atípicos o “outliers”. modo de ejemplo podemos observar el siguiente gráfico de la aportación… donde se aprecia un valor negativo enorme y erróneo sin ninguna duda.
plot(DHI$`APORTACION (m3/s)`,
xlab = "Date",
ylab= "Aportacion [m³/s]")
Investigando sobre la mejor manera de eliminar outliers encontre una función hecha por un fanático. que funciona de la siguiente manera y la verdad que me gusta.
outlierKD <- function(dt, var) {
var_name <- eval(substitute(var),eval(dt))
tot <- sum(!is.na(var_name))
na1 <- sum(is.na(var_name))
m1 <- mean(var_name, na.rm = T)
par(mar = rep(2, 4))
par(mfrow=c(2, 2), oma=c(0,0,3,0))
boxplot(var_name, main="With outliers")
hist(var_name, main="With outliers", xlab=NA, ylab=NA)
outlier <- boxplot.stats(var_name)$out
mo <- mean(outlier)
var_name <- ifelse(var_name %in% outlier, NA, var_name)
boxplot(var_name, main="Without outliers")
hist(var_name, main="Without outliers", xlab=NA, ylab=NA)
title("Outlier Check", outer=TRUE)
na2 <- sum(is.na(var_name))
message("Outliers identified: ", na2 - na1, " from ", tot, " observations")
message("Proportion (%) of outliers: ", (na2 - na1) / tot*100)
message("Mean of the outliers: ", mo)
m2 <- mean(var_name, na.rm = T)
message("Mean without removing outliers: ", m1)
message("Mean if we remove outliers: ", m2)
dt[as.character(substitute(var))] <- invisible(var_name)
assign(as.character(as.list(match.call())$dt), dt, envir = .GlobalEnv)
message("Outliers successfully removed", "\n")
return(invisible(dt))
}
Eliminamos valores negativos y pasamos el filtro para outliers.
DHI$`APORTACION (m3/s)`[which(DHI$`APORTACION (m3/s)`<0)]<- NA
outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 826 from 34552 observations
## Proportion (%) of outliers: 2.39059967585089
## Mean of the outliers: 1075.81016949153
## Mean without removing outliers: 131.482568013429
## Mean if we remove outliers: 108.354577773824
## Outliers successfully removed
outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 317 from 33726 observations
## Proportion (%) of outliers: 0.939927652256419
## Mean of the outliers: 410.450378548896
## Mean without removing outliers: 108.354577773824
## Mean if we remove outliers: 105.488153491574
## Outliers successfully removed
outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 228 from 33409 observations
## Proportion (%) of outliers: 0.682450836600916
## Mean of the outliers: 391.157280701754
## Mean without removing outliers: 105.488153491574
## Mean if we remove outliers: 103.525205991381
## Outliers successfully removed
outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 159 from 33181 observations
## Proportion (%) of outliers: 0.479189897833097
## Mean of the outliers: 380.017798742138
## Mean without removing outliers: 103.525205991381
## Mean if we remove outliers: 102.193901944158
## Outliers successfully removed
outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 168 from 33022 observations
## Proportion (%) of outliers: 0.5087517412634
## Mean of the outliers: 373.650892857143
## Mean without removing outliers: 102.193901944158
## Mean if we remove outliers: 100.805797771961
## Outliers successfully removed
outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 50 from 32854 observations
## Proportion (%) of outliers: 0.152188470201498
## Mean of the outliers: 365.8664
## Mean without removing outliers: 100.805797771961
## Mean if we remove outliers: 100.40179124497
## Outliers successfully removed
outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 4 from 32804 observations
## Proportion (%) of outliers: 0.0121936349225704
## Mean of the outliers: 362.2675
## Mean without removing outliers: 100.40179124497
## Mean if we remove outliers: 100.369856402439
## Outliers successfully removed
outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 1 from 32800 observations
## Proportion (%) of outliers: 0.00304878048780488
## Mean of the outliers: 361.84
## Mean without removing outliers: 100.369856402439
## Mean if we remove outliers: 100.361884508674
## Outliers successfully removed
outlierKD(DHI,`APORTACION (m3/s)`)
## Outliers identified: 0 from 32799 observations
## Proportion (%) of outliers: 0
## Mean of the outliers: NaN
## Mean without removing outliers: 100.361884508674
## Mean if we remove outliers: 100.361884508674
## Outliers successfully removed
plot(DHI$`APORTACION (m3/s)`,
xlab = "Date",
ylab= "Aportacion [m³/s]",
type = "l")
Aplicando esta funcion a los datos de aportación vemos como se han eliminado valores… pero observando los datos… nos damos cuenta de que nuestros datos siguen teniendo demasiado “ruido”. Antes de proseguir hay que rellenar los huecos que hemos ido generando.
Tras hacer la limpieza de los datos, se han generado gran cantidad de NA’s dentro de nuestro dataset.
sum(is.na(DHI$`APORTACION (m3/s)`))
## [1] 7180
A continuación añadimos datos interpolando. Con el paquete imputeTS se puede añadir los datos que faltan interpolando de maneras diferentes:
library(imputeTS)
#Rellenamos Huecos
DHI$`APORTACION (m3/s)`<- na.interpolation(DHI$`APORTACION (m3/s)`)
#Comprobamos
sum(is.na(DHI$`APORTACION (m3/s)`))
## [1] 0
Para eliminarlo aplicamos una moving average. Tenemos en cuenta de que tenemos un total de 96 datos diarios. probamos la Moving Average para varios valores.
library(TTR)
for(n in c(1,seq(48,192,length.out = 4))){
Mavg<- SMA(DHI$`APORTACION (m3/s)`, n)
plot(DHI$DATE, Mavg,
xlab = "Date",
ylab= "Aportacion [m³/s]",
type = "l",
main = paste0("Aportación en Belesar. MA= ", n))
}
Se puede apreciar como de este mode se elimina bastante ruido de la señal de aportación. Pero no sé hasta que punto estamos respetando los datos… Para continuar se decide que utilizaremos una Moving Average con 96 puntos(cantidad de datos generados en 1 dia)… que sería algo parecido a una media diaria (aunque no es lo mismo).
Tras ejecutar la funcion en buscar de outliers vemos que no hay que en este dataset hay menos datos erroneos.
outlierKD(DHI,`NIVEL EMBALSE (msnm)`)
## Outliers identified: 1 from 39978 observations
## Proportion (%) of outliers: 0.00250137575666617
## Mean of the outliers: 0
## Mean without removing outliers: 310.775323928161
## Mean if we remove outliers: 310.783097781224
## Outliers successfully removed
DHI$`NIVEL EMBALSE (msnm)`<- na.interpolation(DHI$`NIVEL EMBALSE (msnm)`)
A continuación comprobaremos la correlacion que existe entre la aportacion y la variacion del nivel de la presa de Belesar.
x<- 96
aportacion_SMA<- SMA(DHI$`APORTACION (m3/s)`,x)
nivel_SMA<- SMA(c(0,diff(DHI$`NIVEL EMBALSE (msnm)`)*8000+150), x)
aportacion_SMA<- aportacion_SMA[!is.na(aportacion_SMA)]
nivel_SMA<- nivel_SMA[!is.na(nivel_SMA)]
plot(aportacion_SMA,type = "l",
xlab =" " ,ylab = "",xaxt= 'n',yaxt='n',
main=paste0("Variacion del nivel y aportacion\n Correlación de: ",
round(cor(aportacion_SMA, nivel_SMA), 3)))
lines(nivel_SMA, col = "red")
Otro análisis previo que se puede hacer es el de la crosscorrelation… para observar con que desfase se obtiene la mejor correlation. Esto es super útil para analizar correlacione de variables que estan desfasadas en el tiempo.
ccf_belesar<- ccf(aportacion_SMA, nivel_SMA, lag.max = 5000)
#Máxima correlación
max(ccf_belesar$acf)
## [1] 0.5892262
#Con cuanto desfase se produce la máxima correlación.
ccf_belesar$lag[which.max(ccf_belesar$acf)]
## [1] 34
De esta manera podemos darnos cuenta de que estas variables obtienen la mejor correlacion (0.58) para un retraso en la aportación de 34 valores, teniendo en cuenta que nuestro dataset es de datos 15 minutales esto supone 9 horas, es decir, los datos de aportación influyen en el nivel de la presa con un retardo de medio día.
A continuación construiremos una tabla con las variables aportacion y nivel obtenidas por diferentes métodos:
#Retrasamos hacia detrás todos los datos de lluvia porque la acumulada de toda la hora la suma en la hora siguiente....
DHI$`LLUVIA ACUMULADA DÍA (l/m2)`<- lead(DHI$`LLUVIA ACUMULADA DÍA (l/m2)`)
Desacumular_lluvia_2018<- DHI[which(year(DHI$DATE)==2018),] %>% group_by(yday(DATE)) %>% mutate(desacumulada= c(diff(`LLUVIA ACUMULADA DÍA (l/m2)`),0),
lluvia=ifelse(desacumulada>=0, desacumulada, 0))
Desacumular_lluvia_2019<- DHI[which(year(DHI$DATE)==2019),] %>% group_by(yday(DATE)) %>% mutate(desacumulada= c(diff(`LLUVIA ACUMULADA DÍA (l/m2)`),0),
lluvia=ifelse(desacumulada>=0, desacumulada, 0))
Desacumular_lluvia<- rbind(Desacumular_lluvia_2018, Desacumular_lluvia_2019)
Medias_horarias_2018<- as.data.frame(Desacumular_lluvia[year(Desacumular_lluvia$DATE)==2018,]) %>% group_by(yday(DATE), hour(DATE)) %>%
summarize(., Acum_horaria=sum(lluvia, na.rm = T),
aport_mean=mean(`APORTACION (m3/s)`, na.rm = T),
nivel_mean=mean(`NIVEL EMBALSE (msnm)`, na.rm = T))
Medias_horarias_2019<- as.data.frame(Desacumular_lluvia[year(Desacumular_lluvia$DATE)==2019,]) %>% group_by(yday(DATE),hour(DATE)) %>%
summarize(., Acum_horaria=sum(lluvia, na.rm = T),
aport_mean=mean(`APORTACION (m3/s)`, na.rm = T),
nivel_mean=mean(`NIVEL EMBALSE (msnm)`, na.rm = T))
Medias_horarias<- rbind(Medias_horarias_2018, Medias_horarias_2019)
vector_Date<- seq(range(DHI$DATE)[1],range(DHI$DATE)[2],
by="hour")
Lluvia_acum_horaria<- as.data.frame(cbind(as.character(vector_Date[2:length(vector_Date)]),
Medias_horarias[,3:5]))
Lluvia_acum_horaria$`as.character(vector_Date[2:length(vector_Date)])`<- ymd_hms(Lluvia_acum_horaria$`as.character(vector_Date[2:length(vector_Date)])`)
colnames(Lluvia_acum_horaria)<- c("Date", "Lluvia_mm", "aport_mean", "nivel_mean")
Comprobacoin de máxima correlacion desplazada para diferentes valores de SMA.
for (n in seq(6,24*3, by=6)) {
y<- n
aportacion_mean_SMA<- SMA(Lluvia_acum_horaria$aport_mean,y)
nivel_mean_SMA<- SMA(Lluvia_acum_horaria$nivel_mean,y)
Tabla_DHI<- as.data.frame(cbind(Lluvia_acum_horaria,
aportacion_mean_SMA,
nivel_mean_SMA))
colnames(Tabla_DHI)<- c(names(Lluvia_acum_horaria), "aport_mean_SMA", "nivel_mean_SMA")
Tabla_DHI<- Tabla_DHI[complete.cases(Tabla_DHI),]
ccf_aport_mean<- ccf(Tabla_DHI$aport_mean_SMA,
Tabla_DHI$nivel_mean_SMA,
lag.max = 5000,
plot = T)
print(paste0("Máxima correlacion de ",
round(max(ccf_aport_mean$acf), digits = 3) ,
" .Para un desfase de: ",
round(ccf_aport_mean$lag[which.max(ccf_aport_mean$acf)]/24, digits = 2),
" días. Con un SMA de :",
n,
" horas"))
}
## [1] "Máxima correlacion de 0.609 .Para un desfase de: -76.79 días. Con un SMA de :6 horas"
## [1] "Máxima correlacion de 0.616 .Para un desfase de: -76.58 días. Con un SMA de :12 horas"
## [1] "Máxima correlacion de 0.621 .Para un desfase de: -76.42 días. Con un SMA de :18 horas"
## [1] "Máxima correlacion de 0.624 .Para un desfase de: -76.38 días. Con un SMA de :24 horas"
## [1] "Máxima correlacion de 0.627 .Para un desfase de: -76.29 días. Con un SMA de :30 horas"
## [1] "Máxima correlacion de 0.629 .Para un desfase de: -76.25 días. Con un SMA de :36 horas"
## [1] "Máxima correlacion de 0.632 .Para un desfase de: -76.21 días. Con un SMA de :42 horas"
## [1] "Máxima correlacion de 0.634 .Para un desfase de: -76.08 días. Con un SMA de :48 horas"
## [1] "Máxima correlacion de 0.636 .Para un desfase de: -75.92 días. Con un SMA de :54 horas"
## [1] "Máxima correlacion de 0.638 .Para un desfase de: -75.88 días. Con un SMA de :60 horas"
## [1] "Máxima correlacion de 0.64 .Para un desfase de: -75.79 días. Con un SMA de :66 horas"
## [1] "Máxima correlacion de 0.642 .Para un desfase de: -75.67 días. Con un SMA de :72 horas"
Se puede observar como casi todos los casos probando diferentes SMA’s llegan a la conclusión de que la correlación máxima se tiene para 3 días y poco entre aportación y nivel… Se logra mejorar la correlacion aumentando el SMA. La mejor correlación aparece para un desfase de entorno a 75 días. Lo cual no tiene mucho sentido. o por lo menos a nosotros no nos vale…
Añadimos ahora los SMA’s de todo el dataset y el SMA sobre la media horaria.
#se selecciona 96 porque es la cantidad de datos que corresponden a 1 día.
x<- 96
aportacion_SMA<- SMA(DHI$`APORTACION (m3/s)`,x)
nivel_SMA<- SMA(DHI$`NIVEL EMBALSE (msnm)`, x)
aportacion_horaria<- aportacion_SMA[DHI$DATE%in%vector_Date]
nivel_horario<- nivel_SMA[DHI$DATE%in%vector_Date]
#Se ponen 24 horas para que equivalgan a 1 día
y<- 24
aportacion_mean_SMA<- SMA(Lluvia_acum_horaria$aport_mean,y)
nivel_mean_SMA<- SMA(Lluvia_acum_horaria$nivel_mean,y)
Tabla_DHI<- as.data.frame(cbind(Lluvia_acum_horaria,
aportacion_horaria[2:length(aportacion_horaria)],
nivel_horario[2:length(nivel_horario)],
aportacion_mean_SMA,
nivel_mean_SMA))
colnames(Tabla_DHI)<- c(names(Lluvia_acum_horaria), "aport_SMA", "nivel_SMA",
"aport_mean_SMA", "nivel_mean_SMA")
Tabla_DHI<- Tabla_DHI[complete.cases(Tabla_DHI),]
A continuación representamos la aportacion y el nivel. para los tres casos
k<- max(Tabla_DHI$nivel_mean_SMA)/max(Tabla_DHI$aport_mean_SMA)
ggplot(data=Tabla_DHI, aes(x=Date))+
geom_line(aes(y=aport_mean_SMA))+
xlab("Date")+ylab("Aportación [m³/s]")+
theme(panel.background = element_blank(),
panel.grid = element_blank())+
geom_line(aes(y=nivel_mean_SMA/k),group=1, color="red") +
scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Nivel [msnm]",
breaks = seq(min(Tabla_DHI$nivel_mean_SMA),
max(Tabla_DHI$nivel_mean_SMA),
by=10)),
breaks = seq(min(Tabla_DHI$aport_mean_SMA),
max(Tabla_DHI$aport_mean_SMA),
by=20)) +
ggtitle("Nivel y aportación. SMA (24 horas) de la media horaria")
ccf_DHI<- ccf(Tabla_DHI$aport_mean_SMA,
Tabla_DHI$nivel_mean_SMA,
lag.max = 1000)
print(paste0("Máxima correlacion de ",
round(max(ccf_DHI$acf), digits = 3) ,
" .Para un desfase de: ",
round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2),
" días"))
## [1] "Máxima correlacion de 0.583 .Para un desfase de: -41.67 días"
k<- max(Tabla_DHI$nivel_SMA)/max(Tabla_DHI$aport_SMA)
ggplot(data=Tabla_DHI, aes(x=Date))+
geom_line(aes(y=aport_SMA))+
xlab("Date")+ylab("Aportación [m³/s]")+
theme(panel.background = element_blank(),
panel.grid = element_blank())+
geom_line(aes(y=nivel_SMA/k),group=1, color="red") +
scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Nivel [msnm]",
breaks = seq(min(Tabla_DHI$nivel_SMA),
max(Tabla_DHI$nivel_SMA),
by=10)),
breaks = seq(min(Tabla_DHI$aport_SMA),
max(Tabla_DHI$aport_SMA),
by=20)) +
ggtitle("Nivel y aportación. SMA (96 datos = 1 día) de los datos 15-minutales")
ccf_DHI<- ccf(Tabla_DHI$aport_SMA,
Tabla_DHI$nivel_SMA,
lag.max = 4000)
print(paste0("Máxima correlacion de ",
round(max(ccf_DHI$acf), digits = 3) ,
" .Para un desfase de: ",
round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2),
" días"))
## [1] "Máxima correlacion de 0.624 .Para un desfase de: -76.38 días"
k<- max(Tabla_DHI$nivel_mean)/max(Tabla_DHI$aport_mean)
ggplot(data=Tabla_DHI, aes(x=Date))+
geom_line(aes(y=aport_mean))+
xlab("Date")+ylab("Aportación [m³/s]")+
theme(panel.background = element_blank(),
panel.grid = element_blank())+
geom_line(aes(y=nivel_mean/k),group=1, color="red") +
scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Nivel [msnm]",
breaks = seq(min(Tabla_DHI$nivel_mean),
max(Tabla_DHI$nivel_mean),
by=10)),
breaks = seq(min(Tabla_DHI$aport_mean),
max(Tabla_DHI$aport_mean),
by=20)) +
ggtitle("Nivel y aportación. Medias horarias")
ccf_DHI<- ccf(Tabla_DHI$aport_mean,
Tabla_DHI$nivel_mean,
lag.max = 1000)
print(paste0("Máxima correlacion de ",
round(max(ccf_DHI$acf), digits = 3) ,
" .Para un desfase de: ",
round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2),
" días"))
## [1] "Máxima correlacion de 0.54 .Para un desfase de: -41.67 días"
Por lo anteriormente expuesto… se llega a la conclusión de que el mejor método para relacionar la aportación y el nivel es haciendo una moving average sobre las medias horarias… un SMA de 142 valores… entorno a 6 dias.
Otra comprobacion que considero pertinente es ver si la lluvia y la aportacion correlan debidamente…
k<- max(Tabla_DHI$Lluvia_mm)/max(Tabla_DHI$aport_mean_SMA)
ggplot(data=Tabla_DHI, aes(x=Date))+
geom_line(aes(y=aport_mean_SMA))+
xlab("Date")+ylab("Aportación [m³/s]")+
theme(panel.background = element_blank(),
panel.grid = element_blank())+
geom_line(aes(y=Lluvia_mm/k),group=1, color="red") +
scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Lluvia Horaria [l/m²]",
breaks = seq(min(Tabla_DHI$Lluvia_mm),
max(Tabla_DHI$Lluvia_mm),
by=10)),
breaks = seq(min(Tabla_DHI$aport_mean_SMA),
max(Tabla_DHI$aport_mean_SMA),
by=20)) +
ggtitle("Lluvia y aportacion")
ccf_DHI<- ccf(Tabla_DHI$aport_mean_SMA,
Tabla_DHI$Lluvia_mm,
lag.max = 1000)
print(paste0("Máxima correlacion de ",
round(max(ccf_DHI$acf), digits = 3) ,
" .Para un desfase de: ",
round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2),
" días"))
## [1] "Máxima correlacion de 0.109 .Para un desfase de: 1.71 días"
Tabla_DHI_cut<- Tabla_DHI[which(month(Tabla_DHI$Date)==12), ]
k<- max(Tabla_DHI_cut$Lluvia_mm)/max(Tabla_DHI_cut$aport_mean_SMA)
ggplot(data=Tabla_DHI_cut, aes(x=Date))+
geom_line(aes(y=aport_mean_SMA))+
xlab("Date")+ylab("Aportación [m³/s]")+
theme(panel.background = element_blank(),
panel.grid = element_blank())+
geom_line(aes(y=Lluvia_mm/k),group=1, color="red") +
scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Nivel [msnm]",
breaks = seq(min(Tabla_DHI_cut$Lluvia_mm),
max(Tabla_DHI_cut$Lluvia_mm),
by=10)),
breaks = seq(min(Tabla_DHI_cut$aport_mean_SMA),
max(Tabla_DHI_cut$aport_mean_SMA),
by=20)) +
ggtitle("Lluvia y aportacion. Diciembre")
ccf_DHI<- ccf(Tabla_DHI_cut$aport_mean_SMA,
Tabla_DHI_cut$Lluvia_mm,
lag.max = 1000)
print(paste0("Máxima correlacion de ",
round(max(ccf_DHI$acf), digits = 3) ,
" .Para un desfase de: ",
round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2),
" días"))
## [1] "Máxima correlacion de 0.371 .Para un desfase de: 1.62 días"
Tabla_DHI_cut<- Tabla_DHI[which(month(Tabla_DHI$Date)==12 | year(Tabla_DHI$Date)==2019), ]
k<- max(Tabla_DHI_cut$Lluvia_mm)/max(Tabla_DHI_cut$aport_mean_SMA)
ggplot(data=Tabla_DHI_cut, aes(x=Date))+
geom_line(aes(y=aport_mean_SMA))+
xlab("Date")+ylab("Aportación [m³/s]")+
theme(panel.background = element_blank(),
panel.grid = element_blank())+
geom_line(aes(y=Lluvia_mm/k),group=1, color="red") +
scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Nivel [msnm]",
breaks = seq(min(Tabla_DHI_cut$Lluvia_mm),
max(Tabla_DHI_cut$Lluvia_mm),
by=10)),
breaks = seq(min(Tabla_DHI_cut$aport_mean_SMA),
max(Tabla_DHI_cut$aport_mean_SMA),
by=20)) +
ggtitle("Lluvia y aportacion. Diciembre-Enero-Febrero")
ccf_DHI<- ccf(Tabla_DHI_cut$aport_mean_SMA,
Tabla_DHI_cut$Lluvia_mm,
lag.max = 1000)
print(paste0("Máxima correlacion de ",
round(max(ccf_DHI$acf), digits = 3) ,
" .Para un desfase de: ",
round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2),
" días"))
## [1] "Máxima correlacion de 0.332 .Para un desfase de: 1.71 días"
Se puede observar como la correlación entre la aportación y la lluvia es apriori muy mala… habrá que hacer trikiñuelas para poder salvar esto… En principio me parece normal que estas variables tengan mala correlación… la aportación del río si que dependerá e la lluvia… pero influyen muchos más factores… además. La lluvia registrada en belesar no tiene porqué ser la mejor para estimar la aportación del río que llega a belesar… puede estar lluviendo muchos kilometros hacia arriba en el valle…
Antes de proseguir guardamos este histórico que hemos construido.
path_dataDHI <- here::here('Data/Parques/Belesar/Historico/DHI_historico_afinado.RDS')
#saveRDS(Tabla_DHI,path_dataDHI)
Los datos del histórico de WRF están organizados en formato lista… hay para cada localización 2 dataframes;
A continuación se presenta el loop y las funciones que se emplean para generar todo el histórico WRF de Belesar… Cabe mencionar que contamos con una carpeta con Archivos RDS ordenados por fechas para la localización de belesar.
Usamos las siguientes funciones.
lon_lat_df_ls<- function(parque_list){
fecha<- names(parque_list)
fecha[!str_detect(fecha, " ")]<- paste0(fecha[!str_detect(fecha, " ")], " 00:00:00")
n_lon<- parque_list[[1]]$lon
n_lat<- parque_list[[1]]$lat
nombres<- paste0(n_lon,"__",n_lat)
list_localizaciones<- list()
for (localizaciones in 1:length(parque_list[[1]][,1]) ) {
l<- lapply(parque_list, function(x) x[localizaciones,])
df <- data.frame(matrix(unlist(l), nrow=length(l), byrow=T))
df<- as.data.frame(cbind(ymd_hms(fecha), df))
colnames(df)<-c("fechas", colnames(parque_list[[1]]))
list_localizaciones[[localizaciones]]<- df
}
names(list_localizaciones)<- nombres
return(list_localizaciones)
}
uv_transformation<- function(tabla_comp){
nombres<- colnames(tabla_comp)
u10<- tabla_comp$U10_MEAN
v10<- tabla_comp$V10_MEAN
u10_max<- tabla_comp$U10_MAX
v10_max<- tabla_comp$V10_MAX
wind_abs <- sqrt(u10^2 + v10^2)
wind_dir_rad <- atan2(u10/wind_abs, v10/wind_abs)
wind_dir_deg1 <- wind_dir_rad * 180/pi
wind_dir_deg2 <- wind_dir_deg1+ 180
wind_abs_max <- sqrt(u10_max^2 + v10_max^2)
wind_dir_rad_max <- atan2(u10_max/wind_abs_max, v10_max/wind_abs_max)
wind_dir_deg1_max <- wind_dir_rad_max * 180/pi
wind_dir_deg2_max <- wind_dir_deg1_max+ 180
tabla_comp<- as.data.frame(cbind(tabla_comp,wind_abs,wind_dir_deg2,
wind_abs_max,wind_dir_deg2_max))
colnames(tabla_comp)<- c(nombres,"WS","WD","WS_MAX","WD_MAX")
tabla_comp$U10_MAX<- NULL
tabla_comp$V10_MAX<- NULL
tabla_comp$U10_MEAN<- NULL
tabla_comp$V10_MEAN<- NULL
return(tabla_comp)
}
extract_rain_data<- function(Belesar_lolat_df){
rain<- Belesar_lolat_df[,c(1,2,3,4,5)]
rain$pre_acum<- rain$RAINC+rain$RAINNC
rain$RAINC<- NULL
rain$RAINNC<- NULL
prep_hourly<- vector()
for (i in 1:length(rain$pre_acum)) {
if(i==1){prep_hourly[i]<- rain$pre_acum[i]}else{
prep_hourly[i]<- rain$pre_acum[i]-rain$pre_acum[i-1]
}
}
rain$prep_hourly<- prep_hourly
return(rain)
}
extract_rain_data2<- function(Belesar_lolat_df){
rain<- Belesar_lolat_df[,c(1,2,3,4,5,6,10,17,20)]
rain$pre_acum<- rain$RAINC+rain$RAINNC
prep_hourly<- vector()
for (i in 1:length(rain$pre_acum)) {
if(i==1){prep_hourly[i]<- rain$pre_acum[i]}else{
prep_hourly[i]<- rain$pre_acum[i]-rain$pre_acum[i-1]
}
}
rain$prep_hourly<- prep_hourly
return(rain)
}
Return_periodo_Belesar<- function(){
All_files_Belesar<- list.files(here::here('Data/Parques/Belesar/'),
recursive = F, full.names = T)
RDS_Belesar<- All_files_Belesar[str_detect(All_files_Belesar, ".RDS")]
RDS_Belesar1<- RDS_Belesar[!str_detect(RDS_Belesar, "E001")]
Belesar_data<- readRDS(RDS_Belesar1[1])
Belesar_lolat<- lon_lat_df_ls(Belesar_data)
Belesar_lolat1<- lapply(Belesar_lolat, uv_transformation)
Belesar_rain<- lapply(Belesar_lolat1, extract_rain_data)
fecha_ini<- Belesar_rain$`-8.02328491210938__42.1343421936035`$fechas[1]
Belesar_data<- readRDS(RDS_Belesar1[length(RDS_Belesar)])
Belesar_lolat<- lon_lat_df_ls(Belesar_data)
Belesar_lolat1<- lapply(Belesar_lolat, uv_transformation)
Belesar_rain<- lapply(Belesar_lolat1, extract_rain_data)
fecha_last<- Belesar_rain$`-8.02328491210938__42.1343421936035`$fechas[length(Belesar_rain$`-8.02328491210938__42.1343421936035`$fechas)]
periodo_WRF<- seq(fecha_ini, fecha_last, by="hour")
Tabla_WRF<- as.data.frame(matrix(ncol = 10, nrow = length(periodo_WRF)))
colnames(Tabla_WRF)<- colnames(Belesar_rain[[1]])
Tabla_WRF$fechas<- periodo_WRF
return(Tabla_WRF)
}
#Listamos archivos dentro de la carpeta de Belesar
All_files_Belesar<- list.files(here::here('Data/Parques/Belesar/'),
recursive = F, full.names = T)
#Detectamos cuales son RDS
RDS_Belesar<- All_files_Belesar[str_detect(All_files_Belesar, ".RDS")]
#Eliminamos los RDS que no son de WRF
RDS_Belesar1<- RDS_Belesar[!str_detect(RDS_Belesar, "E001")]
#Construimos Lista para cada instante de tiempo
Lista_localizacion<- list()
for (i in 1:length(RDS_Belesar1)) {
Belesar_data<- readRDS(RDS_Belesar1[i])
Belesar_lolat<- lon_lat_df_ls(Belesar_data)
Belesar_lolat1<- lapply(Belesar_lolat, uv_transformation)
Belesar_rain<- lapply(Belesar_lolat1, extract_rain_data2)
Lista_localizacion[[i]]<- Belesar_rain
}
#Nombramos la lista
names_fechas<- sapply(RDS_Belesar1, function(x){
r<- str_split(x, "/")
return(str_remove(str_remove(r[[1]][length(r[[1]])], ".RDS"), "Belesar_"))
})
names(Lista_localizacion)<- names_fechas
#Creamos una lista por localización
Lista_localizacion2<- list()
for (i in 1:length(Lista_localizacion[[1]])) {
Lista_localizacion2[[i]]<- lapply(Lista_localizacion,
function(x) return(x[[i]]))
}
#Nombramos la lista
names(Lista_localizacion2)<- names(Lista_localizacion[[1]])
#Guardamos la lista
path_lista_total<- here::here('Data/Parques/Belesar/Historico/')
nombre_lista<- paste0(path_lista_total, 'Historico_WRF_Belesar_Variables.RDS')
#saveRDS(Lista_localizacion2, nombre_lista)
#Cargamos lista
path_lista_total<- here::here('Data/Parques/Belesar/Historico/')
nombre_lista<- paste0(path_lista_total, 'Historico_WRF_Belesar_Variables.RDS')
Lista_total1<- readRDS(nombre_lista)
#Juntamos todos los Dataframes
Lista_total_MF<- lapply(Lista_total1, function(x) bind_rows(x))
#creamos dos data.frames... uno para D1 y otro para D2
Lista_d1_d2_loc<- list()
for (i in 1:length(Lista_total_MF)) {
p<- Lista_total_MF[[i]]
d1<- p[duplicated(p$fechas),]
d1<-d1[!duplicated(d1$fechas),]
d2<- p[!duplicated(p$fechas),]
d2_qneed1<-d2[!(d2$fechas%in%d1$fechas),]
d1_2<-bind_rows(d1,d2_qneed1)
d1_2<-d1_2[order(d1_2$fechas),]
d2<-d2[order(d2$fechas),]
d1_2$pre_acum<- NULL
d2$pre_acum<- NULL
colnames(d1_2)<- c("Date", "LON", "LAT", "RAINC", "RAINNC","RAINSH", "T02_MEAN","PSFC","WS_MAX", "prep_hourly")
colnames(d2)<- c("Date", "LON", "LAT", "RAINC", "RAINNC","RAINSH", "T02_MEAN","PSFC","WS_MAX", "prep_hourly")
lista_loc_d12<- list(d1_2,d2)
names(lista_loc_d12)<- c("D1", "D2")
Lista_d1_d2_loc[[i]]<- lista_loc_d12
}
#nombramos la lsta
names(Lista_d1_d2_loc)<- names(Lista_total_MF)
Tabla_periodo1<- Return_periodo_Belesar()
colnames(Tabla_periodo1)<- c("Date", "LON", "LAT", "RAINC", "RAINNC","RAINSH", "T02_MEAN","PSFC","WS_MAX", "prep_hourly")
Lista_d1_d2_loc2<- list()
for (j in 1:length(Lista_d1_d2_loc)) {
prueba_list<- Lista_d1_d2_loc[[j]]
lista_retorno<- list()
for(i in 1:2){
prueba<- prueba_list[[i]]
Tabla_periodo<- Tabla_periodo1
Tabla_periodo$LON[match(prueba$Date,Tabla_periodo$Date)] <- prueba$LON
Tabla_periodo$LAT[match(prueba$Date,Tabla_periodo$Date)] <- prueba$LAT
Tabla_periodo$RAINC[match(prueba$Date,Tabla_periodo$Date)] <- prueba$RAINC
Tabla_periodo$RAINNC[match(prueba$Date,Tabla_periodo$Date)] <- prueba$RAINNC
Tabla_periodo$RAINSH[match(prueba$Date,Tabla_periodo$Date)] <- prueba$RAINSH
Tabla_periodo$T02_MEAN[match(prueba$Date,Tabla_periodo$Date)] <- prueba$T02_MEAN
Tabla_periodo$PSFC[match(prueba$Date,Tabla_periodo$Date)] <- prueba$PSFC
Tabla_periodo$WS_MAX[match(prueba$Date,Tabla_periodo$Date)] <- prueba$WS_MAX
Tabla_periodo$prep_hourly[match(prueba$Date,Tabla_periodo$Date)] <- prueba$prep_hourly
lista_retorno[[i]]<- Tabla_periodo
}
names(lista_retorno)<- c("D1", "D2")
Lista_d1_d2_loc2[[j]]<- lista_retorno
}
names(Lista_d1_d2_loc2)<- names(Lista_d1_d2_loc)
#Creamos lista con las variables afinadas
Lista_d1_d2_loc3<- lapply(Lista_d1_d2_loc2, function(x){
x[[1]]$RAINC<- NULL
x[[1]]$RAINNC<- NULL
x[[1]]$RAINSH<- NULL
x[[2]]$RAINC<- NULL
x[[2]]$RAINNC<- NULL
x[[2]]$RAINSH<- NULL
return(x)})
names(Lista_d1_d2_loc3)<- names(Lista_d1_d2_loc2)
#Guardamos
path_hist_WRF<- here::here('Data/Parques/Belesar/Historico/Historico_WRF_Belesar_Variables_D1D2.RDS')
#saveRDS(Lista_d1_d2_loc3,path_hist_WRF)
Juntamos a cada data.frame de WRF la tabla afinada de DHI.
#Usamos dplyr para juntar ambos datasets
Belesar_WRF<- readRDS(here::here('Data/Parques/Belesar/Historico/Historico_WRF_Belesar_Variables_D1D2.RDS'))
Belesar_DHI<- readRDS(here::here('Data/Parques/Belesar/Historico/DHI_historico_afinado.RDS'))
df2<- Belesar_DHI
Belesar_Merge<- list()
for (j in 1:length(Belesar_WRF)) {
lista_retorno<- list()
for(i in 1:2){
df1<- Belesar_WRF[[j]][[i]]
Merge_table<- left_join(df1, df2, by=c("Date"))
lista_retorno[[i]]<- Merge_table
}
names(lista_retorno)<- c("D1", "D2")
Belesar_Merge[[j]]<- lista_retorno
}
names(Belesar_Merge)<- names(Belesar_WRF)
#Belesar merge completecases
Belesar_Merge_cc<- list()
for (j in 1:length(Belesar_Merge)) {
lista_retorno<- list()
for(i in 1:2){
df1<- Belesar_Merge[[j]][[i]]
df1$Acumulated<- NULL
Table_fine<- df1[complete.cases(df1),]
lista_retorno[[i]]<- Table_fine
}
names(lista_retorno)<- c("D1", "D2")
Belesar_Merge_cc[[j]]<- lista_retorno
}
names(Belesar_Merge_cc)<- names(Belesar_Merge)
#Guardamos
path_hist_WRF<- here::here('Data/Parques/Belesar/Historico/Hist_D1D2_DHI_MERGED.RDS')
#saveRDS(Belesar_Merge_cc,path_hist_WRF)
#importar
path_hist_WRF<- here::here('Data/Parques/Belesar/Historico/Hist_D1D2_DHI_MERGED.RDS')
clean_data<- readRDS(path_hist_WRF)
#cortar en entrenamiento y predicción
cut_train<- lapply(clean_data, function(x){
y<- lapply(x, function(r){
fecha_ini<- ymd("2018/10/01")
fecha_end<- ymd("2019/02/01")
Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
return(Jan_data)
})
return(y)
})
cut_predict<- lapply(clean_data, function(x){
y<- lapply(x, function(r){
fecha_ini<- ymd("2019/02/01")
fecha_end<- ymd("2019/02/20")
Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
return(Jan_data)
})
return(y)
})
lista_TL<- list()
for (i in 1:length(cut_train)) {
lista_d1d2<- list()
for (j in 1:2) {
data_predict<- cut_train[[i]][[j]]
data_predict2<- cut_predict[[i]][[j]]
fit_1 <- lm(Lluvia_mm ~ prep_hourly, data = data_predict)
fit_2 <- lm(Lluvia_mm ~ prep_hourly + T02_MEAN, data = data_predict)
fit_3 <- lm(Lluvia_mm ~ prep_hourly + PSFC, data = data_predict)
fit_4 <- lm(Lluvia_mm ~ prep_hourly + WS_MAX, data = data_predict)
fit_5 <- lm(Lluvia_mm ~ prep_hourly * T02_MEAN, data = data_predict)
fit_6 <- lm(Lluvia_mm ~ prep_hourly * PSFC, data = data_predict)
fit_7 <- lm(Lluvia_mm ~ prep_hourly * WS_MAX, data = data_predict)
uncorrected<-data_predict2$prep_hourly
prediction_rain<- predict(fit_1, data.frame(prep_hourly =data_predict2$prep_hourly))
prediction_rain2<- predict(fit_2, data.frame(prep_hourly =data_predict2$prep_hourly,
T02_MEAN =data_predict2$T02_MEAN))
prediction_rain3<- predict(fit_3, data.frame(prep_hourly =data_predict2$prep_hourly,
PSFC =data_predict2$PSFC))
prediction_rain4<- predict(fit_4, data.frame(prep_hourly =data_predict2$prep_hourly,
WS_MAX =data_predict2$WS_MAX))
prediction_rain5<- predict(fit_5, data.frame(prep_hourly =data_predict2$prep_hourly,
T02_MEAN =data_predict2$T02_MEAN))
prediction_rain6<- predict(fit_6, data.frame(prep_hourly =data_predict2$prep_hourly,
PSFC =data_predict2$PSFC))
prediction_rain7<- predict(fit_7, data.frame(prep_hourly =data_predict2$prep_hourly,
WS_MAX =data_predict2$WS_MAX))
observed_rain<-data_predict2$Lluvia_mm
tabla_cor<-cbind(cor(uncorrected, observed_rain),
cor(prediction_rain, observed_rain),
cor(prediction_rain2, observed_rain),
cor(prediction_rain3, observed_rain),
cor(prediction_rain4, observed_rain),
cor(prediction_rain5, observed_rain),
cor(prediction_rain6, observed_rain),
cor(prediction_rain7, observed_rain))
fit_1 <- svm(Lluvia_mm ~ prep_hourly, data = data_predict)
fit_2 <- svm(Lluvia_mm ~ prep_hourly + T02_MEAN, data = data_predict)
fit_3 <- svm(Lluvia_mm ~ prep_hourly + PSFC, data = data_predict)
fit_4 <- svm(Lluvia_mm ~ prep_hourly + WS_MAX, data = data_predict)
fit_5 <- svm(Lluvia_mm ~ prep_hourly * T02_MEAN, data = data_predict)
fit_6 <- svm(Lluvia_mm ~ prep_hourly * PSFC, data = data_predict)
fit_7 <- svm(Lluvia_mm ~ prep_hourly * WS_MAX, data = data_predict)
uncorrected<-data_predict2$prep_hourly
prediction_rain<- predict(fit_1, data.frame(prep_hourly =data_predict2$prep_hourly))
prediction_rain2<- predict(fit_2, data.frame(prep_hourly =data_predict2$prep_hourly,
T02_MEAN =data_predict2$T02_MEAN))
prediction_rain3<- predict(fit_3, data.frame(prep_hourly =data_predict2$prep_hourly,
PSFC =data_predict2$PSFC))
prediction_rain4<- predict(fit_4, data.frame(prep_hourly =data_predict2$prep_hourly,
WS_MAX =data_predict2$WS_MAX))
prediction_rain5<- predict(fit_5, data.frame(prep_hourly =data_predict2$prep_hourly,
T02_MEAN =data_predict2$T02_MEAN))
prediction_rain6<- predict(fit_6, data.frame(prep_hourly =data_predict2$prep_hourly,
PSFC =data_predict2$PSFC))
prediction_rain7<- predict(fit_7, data.frame(prep_hourly =data_predict2$prep_hourly,
WS_MAX =data_predict2$WS_MAX))
observed_rain<-data_predict2$Lluvia_mm
tabla_cor<-cbind(tabla_cor,
cor(prediction_rain, observed_rain),
cor(prediction_rain2, observed_rain),
cor(prediction_rain3, observed_rain),
cor(prediction_rain4, observed_rain),
cor(prediction_rain5, observed_rain),
cor(prediction_rain6, observed_rain),
cor(prediction_rain7, observed_rain))
colnames(tabla_cor)<- c("uncorrected", "1","2", "3","4","5","6","7",
"svm1","svm2", "svm3","svm4","svm5","svm6","svm7")
lista_d1d2[[j]]<- tabla_cor
}
names(lista_d1d2)<- c("D1", "D2")
lista_TL[[i]]<- lista_d1d2
}
colnames(tabla_cor)<- c("uncorrected", "1","2", "3","4","5","6","7",
"svm1","svm2", "svm3","svm4","svm5","svm6","svm7")
names(lista_d1d2)<- c("D1", "D2")
Lista_TL_names<-lapply(lista_TL, function(x){
y<- lapply(x, function(r) {
r<- as.data.frame(r)
names(r)<- c("uncorrected", "1","2", "3","4","5","6","7",
"svm1","svm2", "svm3","svm4","svm5","svm6","svm7")
return(r)
})
names(y)<- c("D1", "D2")
return(y)
})
names(Lista_TL_names)<- names(cut_predict)
x<- lapply(Lista_TL_names, function(x){bind_rows(as.data.frame(x))})
y<- bind_rows(x, .id="id")
cor_table<- as.data.frame(matrix(ncol = 3, nrow = 43))
colnames(cor_table)<- c("id","corr", "name")
for (i in 2:length(y[1,])) {
f<- as.data.frame(cbind(y[which.max(y[,i]), c(1,i)], names( y[which.max(y[,i]), c(1,i)])[2]))
colnames(f)<- c("id","corr", "name")
cor_table<- rbind(cor_table,f )
}
cor_table<- cor_table[complete.cases(cor_table),]
punto_Belesar<- c(42.628577,-7.713948)
extract_lat_lon<- lapply(str_split(cor_table$id,"__"), function(x){
return(as.data.frame(cbind(x[[1]],x[[2]])))
})
extract_lat_lon<- as.data.frame(bind_rows(extract_lat_lon))
cor_table<- as.data.frame(cbind(extract_lat_lon, cor_table[,2:3]))
colnames(cor_table)<- c("LON","LAT", "Corr", "Method")
#Librari Geosphere para medir la distancia entre dos puntos a partir de sus coordenadas
library(geosphere)
vec_dist<- vector()
for (i in 1:length(cor_table[,1])) {
vec_dist[i]<- distm(c(-7.713948, 42.628577), c(as.numeric(cor_table[i,1]),
as.numeric(cor_table[i,2])),
fun = distHaversine)
}
vec_dist<- round(vec_dist/1000, digits = 1)
cor_table<- as.data.frame(cbind(cor_table, vec_dist))
cor_table$LON<-round(as.numeric(cor_table$LON), digits = 2)
cor_table$LAT<-round(as.numeric(cor_table$LAT), digits = 2)
cor_table$Corr<-round(as.numeric(cor_table$Corr), digits = 3)
cor_table
De esta manera llegamos a la conclusión de que empleando diferentes variables entregadas por WRF, tales como presión superficial, viento, y temperatura. No existe una regresión lineal o SVM que pueda mejorar la predicción de la precipitación ofrecida por WRF.
Pasamos ahora a analizar la lluvia ofrecida por WRF comparada con la registrada en Belesar. Ya hemos visto que la correlación mayor se obtiene para un punto a 23 km de distancia.
#Empleamos los datos a partir de octubre y eliminamos datos negativos de
# precipitacion WRF que hemos visto que hay alguno.
clean_data_oct<-lapply(clean_data, function(x){
y<- lapply(x, function(r) {
r<- r[which(r$Date > ymd("2018/10/01")),]
r$prep_hourly<- ifelse(r$prep_hourly < 0,0,r$prep_hourly)
return(r)
})
})
Cor_rain_place<- data.frame(matrix(ncol = 4))
colnames(Cor_rain_place)<- c("LON", "LAT", "Corr", "Dist")
for (i in 1:length(clean_data_oct)) {
Corr<- cor(clean_data_oct[[i]][[1]]$prep_hourly,
clean_data_oct[[i]][[1]]$Lluvia_mm)
LON<- as.numeric(str_split(names(clean_data_oct)[i], "_")[[1]][1])
LAT<- as.numeric(str_split(names(clean_data_oct)[i], "_")[[1]][3])
Dist<- distm(c(-7.713948, 42.628577),
c(as.numeric(LON),
as.numeric(LAT)),
fun = distHaversine)
Cor_rain_place[i,]<- as.data.frame(cbind(LON, LAT, Corr, Dist))
}
Cor_order<- Cor_rain_place[order(Cor_rain_place$Corr, decreasing = T),]
Cor_order
#Plot_rain3 necesita columna Date y columna rain
plot_rain3<- function(data,data2, titulo){
ggplot(data=data, aes(x=Date))+
geom_line(aes(y=rain), stat="identity")+
xlab("Date")+ylab("Lluvia por hora [mm/h]")+theme(panel.background = element_blank(),
panel.grid = element_blank()) +
geom_line(data= data2, aes(x=Date, y=rain), color="red", alpha=0.5)+
geom_text(aes(as.POSIXct(ymd("2019/02/01")),10, label=paste0("Cor: ",round(cor(data$rain,data2$rain),
digits = 2)))) +
ggtitle(titulo) +
theme(plot.title = element_text(hjust = 0.5))
}
#Ploteamos los 5 puntos con mejor correlacion
for (i in order(Cor_rain_place$Corr, decreasing = T)[1:5]) {
data1<- as.data.frame(cbind(as.character(clean_data_oct[[i]][[1]]$Date),
clean_data_oct[[i]][[1]]$prep_hourly))
names(data1)<- c("Date", "rain")
data1$Date<- ymd_hms(data1$Date)
data1$rain<- as.numeric(as.character(data1$rain))
data2<- as.data.frame(cbind(as.character(clean_data_oct[[i]][[1]]$Date),
clean_data_oct[[i]][[1]]$Lluvia_mm))
names(data2)<- c("Date", "rain")
data2$Date<- ymd_hms(data2$Date)
data2$rain<- as.numeric(as.character(data2$rain))
x<- plot_rain3(data = data1 ,
data2 = data2,
titulo = paste0(paste(round(as.numeric(str_split(names(clean_data_oct)[i],
"_")[[1]][c(1,3)]),
digits = 3),
collapse = " "), "\n n = ",
i, " Dist = ",
round((Cor_rain_place$Dist[i]/1000), digits = 2), " km" ))
print(x)
}
Es reseñable que disminuyendo el tiempo de ploteo la correlación aumenta.
clean_data_jan<-lapply(clean_data_oct, function(x){
y<- lapply(x, function(r) r[which(r$Date > ymd("2019/01/01")),])
})
Cor_rain_place<- data.frame(matrix(ncol = 4))
colnames(Cor_rain_place)<- c("LON", "LAT", "Corr", "Dist")
for (i in 1:length(clean_data_jan)) {
Corr<- cor(clean_data_jan[[i]][[1]]$prep_hourly,
clean_data_jan[[i]][[1]]$Lluvia_mm)
LON<- as.numeric(str_split(names(clean_data_jan)[i], "_")[[1]][1])
LAT<- as.numeric(str_split(names(clean_data_jan)[i], "_")[[1]][3])
Dist<- distm(c(-7.713948, 42.628577),
c(as.numeric(LON),
as.numeric(LAT)),
fun = distHaversine)
Cor_rain_place[i,]<- as.data.frame(cbind(LON, LAT, Corr, Dist))
}
Cor_order<- Cor_rain_place[order(Cor_rain_place$Corr, decreasing = T),]
Cor_order<- Cor_rain_place[order(Cor_rain_place$Corr, decreasing = T),]
Cor_order
#Ploteamos los 5 puntos con mejor correlacion
for (i in order(Cor_rain_place$Corr, decreasing = T)[1:5]) {
data1<- as.data.frame(cbind(as.character(clean_data_jan[[i]][[1]]$Date),
clean_data_jan[[i]][[1]]$prep_hourly))
names(data1)<- c("Date", "rain")
data1$Date<- ymd_hms(data1$Date)
data1$rain<- as.numeric(as.character(data1$rain))
data2<- as.data.frame(cbind(as.character(clean_data_jan[[i]][[1]]$Date),
clean_data_jan[[i]][[1]]$Lluvia_mm))
names(data2)<- c("Date", "rain")
data2$Date<- ymd_hms(data2$Date)
data2$rain<- as.numeric(as.character(data2$rain))
x<- plot_rain3(data = data1 ,
data2 = data2,
titulo = paste0(paste(round(as.numeric(str_split(names(clean_data_jan)[i],
"_")[[1]][c(1,3)]),
digits = 3),
collapse = " "), "\n n = ",
i, " Dist = ",
round((Cor_rain_place$Dist[i]/1000), digits = 2), " km" ))
print(x)
}
Apuramos más el rango para que únicamente coja un periodo lluvioso como fue el de el 15 de enero al pocos de febrero
clean_data_jan<-lapply(clean_data_oct, function(x){
y<- lapply(x, function(r) r[which(r$Date > ymd("2019/01/15") & r$Date < ymd("2019/02/10")),])
})
Cor_rain_place<- data.frame(matrix(ncol = 4))
colnames(Cor_rain_place)<- c("LON", "LAT", "Corr", "Dist")
for (i in 1:length(clean_data_jan)) {
Corr<- cor(clean_data_jan[[i]][[1]]$prep_hourly,
clean_data_jan[[i]][[1]]$Lluvia_mm)
LON<- as.numeric(str_split(names(clean_data_jan)[i], "_")[[1]][1])
LAT<- as.numeric(str_split(names(clean_data_jan)[i], "_")[[1]][3])
Dist<- distm(c(-7.713948, 42.628577),
c(as.numeric(LON),
as.numeric(LAT)),
fun = distHaversine)
Cor_rain_place[i,]<- as.data.frame(cbind(LON, LAT, Corr, Dist))
}
Cor_order<- Cor_rain_place[order(Cor_rain_place$Corr, decreasing = T),]
Cor_order
#Ploteamos los 5 puntos con mejor correlacion
for (i in order(Cor_rain_place$Corr, decreasing = T)[1:5]) {
data1<- as.data.frame(cbind(as.character(clean_data_jan[[i]][[1]]$Date),
clean_data_jan[[i]][[1]]$prep_hourly))
names(data1)<- c("Date", "rain")
data1$Date<- ymd_hms(data1$Date)
data1$rain<- as.numeric(as.character(data1$rain))
data2<- as.data.frame(cbind(as.character(clean_data_jan[[i]][[1]]$Date),
clean_data_jan[[i]][[1]]$Lluvia_mm))
names(data2)<- c("Date", "rain")
data2$Date<- ymd_hms(data2$Date)
data2$rain<- as.numeric(as.character(data2$rain))
x<- plot_rain3(data = data1 ,
data2 = data2,
titulo = paste0(paste(round(as.numeric(str_split(names(clean_data_jan)[i],
"_")[[1]][c(1,3)]),
digits = 3),
collapse = " "), "\n n = ",
i, " Dist = ",
round((Cor_rain_place$Dist[i]/1000), digits = 2), " km" ))
print(x)
}
A continuacion se realizan las medias diarias de todos los datos y se comprueba la correlación entre los datos de lluvia y aportación.
Daiyli_avgs<- lapply(clean_data_oct, function(x){
x[[1]] %>% group_by(year(Date), yday(Date)) %>%
summarize(rWRF_acum= sum(prep_hourly),
rDHI_acum= sum(Lluvia_mm),
rWRF_mean= sum(prep_hourly),
rDHI_mean= sum(Lluvia_mm),
aport_SMA_diario=mean(aport_SMA),
aport_mean_diario=mean(aport_mean),
aport_mean_SMA_diario=mean(aport_mean_SMA),
nivel_SMA_diario=mean(nivel_SMA),
nivel_mean_diario=mean(nivel_mean),
nivel_mean_SMA_diario=mean(nivel_mean_SMA))
})
Daiyli_avgs2<- lapply(Daiyli_avgs, function(x){
x$Date<- ymd(ifelse(x$`year(Date)`==2018,
as.character(as.Date(x$`yday(Date)`, origin= "2018-01-01")),
as.character(as.Date(x$`yday(Date)`, origin= "2019-01-01"))))
x$`year(Date)`<- NULL
x$`yday(Date)`<- NULL
return(x)
})
Cor_rain_place<- data.frame(matrix(ncol = 11))
for (i in 1:length(Daiyli_avgs2)) {
Corr1<- cor(Daiyli_avgs2[[i]]$rWRF_acum,
Daiyli_avgs2[[i]]$aport_SMA_diario)
Corr2<- cor(Daiyli_avgs2[[i]]$rWRF_acum,
Daiyli_avgs2[[i]]$aport_mean_diario)
Corr3<- cor(Daiyli_avgs2[[i]]$rWRF_acum,
Daiyli_avgs2[[i]]$aport_mean_SMA_diario)
Corr4<- cor(Daiyli_avgs2[[i]]$rWRF_acum,
Daiyli_avgs2[[i]]$rDHI_acum)
Corr5<- cor(Daiyli_avgs2[[i]]$rWRF_mean,
Daiyli_avgs2[[i]]$aport_SMA_diario)
Corr6<- cor(Daiyli_avgs2[[i]]$rWRF_mean,
Daiyli_avgs2[[i]]$aport_mean_diario)
Corr7<- cor(Daiyli_avgs2[[i]]$rWRF_mean,
Daiyli_avgs2[[i]]$aport_mean_SMA_diario)
Corr8<- cor(Daiyli_avgs2[[i]]$rWRF_mean,
Daiyli_avgs2[[i]]$rDHI_mean)
LON<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][1])
LAT<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][3])
Dist<- distm(c(-7.713948, 42.628577),
c(as.numeric(LON),
as.numeric(LAT)),
fun = distHaversine)
Cor_rain_place[i,]<- as.data.frame(cbind(LON, LAT, round(Corr1, digits = 2),
round(Corr2, digits = 2),
round(Corr3, digits = 2),
round(Corr4, digits = 2),
round(Corr5, digits = 2),
round(Corr6, digits = 2),
round(Corr7, digits = 2),
round(Corr, digits = 2),
Dist))
}
colnames(Cor_rain_place)<- c("LON","LAT","AWRF_aport1",
"AWRF_aport2", "AWRF_aport3",
"AWRF_ADHI","MWRF_aport1",
"MWRF_aport2", "MWRF_aport3",
"MWRF_MDHI","Dist")
Se aprecia que los datos de aportación son muy malos
Cor_order<- Cor_rain_place[order(Cor_rain_place$Dist),]
Cor_order
Cor_rain_place_ccf<- data.frame(matrix(ncol = 11))
for (i in 1:length(Daiyli_avgs2)) {
ccfr1<- ccf(Daiyli_avgs2[[i]]$rWRF_acum,
Daiyli_avgs2[[i]]$aport_SMA_diario, plot = FALSE)
ccfr2<- ccf(Daiyli_avgs2[[i]]$rWRF_acum,
Daiyli_avgs2[[i]]$aport_mean_diario, plot = FALSE)
ccfr3<- ccf(Daiyli_avgs2[[i]]$rWRF_acum,
Daiyli_avgs2[[i]]$aport_mean_SMA_diario, plot = FALSE)
ccfr4<- ccf(Daiyli_avgs2[[i]]$rWRF_acum,
Daiyli_avgs2[[i]]$rDHI_acum, plot = FALSE)
ccfr5<- ccf(Daiyli_avgs2[[i]]$rWRF_mean,
Daiyli_avgs2[[i]]$aport_SMA_diario, plot = FALSE)
ccfr6<- ccf(Daiyli_avgs2[[i]]$rWRF_mean,
Daiyli_avgs2[[i]]$aport_mean_diario, plot = FALSE)
ccfr7<- ccf(Daiyli_avgs2[[i]]$rWRF_mean,
Daiyli_avgs2[[i]]$aport_mean_SMA_diario, plot = FALSE)
ccfr8<- ccf(Daiyli_avgs2[[i]]$rWRF_mean,
Daiyli_avgs2[[i]]$rDHI_mean, plot = FALSE)
LON<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][1])
LAT<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][3])
Dist<- distm(c(-7.713948, 42.628577),
c(as.numeric(LON),
as.numeric(LAT)),
fun = distHaversine)
Cor_rain_place_ccf[i,]<- as.data.frame(cbind(LON, LAT,
round(max(ccfr1$acf), digits = 2),
round(max(ccfr2$acf), digits = 2),
round(max(ccfr3$acf), digits = 2),
round(max(ccfr4$acf), digits = 2),
round(max(ccfr5$acf), digits = 2),
round(max(ccfr6$acf), digits = 2),
round(max(ccfr7$acf), digits = 2),
round(max(ccfr8$acf), digits = 2), Dist))
print(paste("Mejor corr para un retardo de: ",
ccfr1$lag[which.max(ccfr1$acf)],
ccfr2$lag[which.max(ccfr2$acf)],
ccfr3$lag[which.max(ccfr3$acf)],
ccfr4$lag[which.max(ccfr4$acf)],
ccfr5$lag[which.max(ccfr5$acf)],
ccfr6$lag[which.max(ccfr6$acf)],
ccfr7$lag[which.max(ccfr7$acf)],
ccfr8$lag[which.max(ccfr8$acf)], sep = " "))
}
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -6 -5 -6 0 -6 -5 -6 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -3 -4 -3 0 -3 -4 -3 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -2 0 -3 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -2 0 -3 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -13 -12 -13 0 -13 -12 -13 0"
colnames(Cor_rain_place_ccf)<- c("LON","LAT","AWRF_aport1",
"AWRF_aport2", "AWRF_aport3",
"AWRF_ADHI","MWRF_aport1",
"MWRF_aport2","MWRF_aport3",
"MWRF_MDHI","Dist")
Cor_order<- Cor_rain_place_ccf[order(Cor_rain_place_ccf$Dist),]
Cor_order
Antes estabamos usando los datos desde octubre hasta febrero… Ahora empleamos solo los datos desde diciembre hasta febrero.
clean_data_dec<-lapply(clean_data, function(x){
y<- lapply(x, function(r) r[which(r$Date > ymd("2018/12/01")),])
})
Daiyli_avgs<- lapply(clean_data_dec, function(x){
x[[1]] %>% group_by(year(Date), yday(Date)) %>%
summarize(rWRF_acum= sum(prep_hourly),
rDHI_acum= sum(Lluvia_mm),
rWRF_mean= sum(prep_hourly),
rDHI_mean= sum(Lluvia_mm),
aport_SMA_diario=mean(aport_SMA),
aport_mean_diario=mean(aport_mean),
aport_mean_SMA_diario=mean(aport_mean_SMA),
nivel_SMA_diario=mean(nivel_SMA),
nivel_mean_diario=mean(nivel_mean),
nivel_mean_SMA_diario=mean(nivel_mean_SMA))
})
Daiyli_avgs2<- lapply(Daiyli_avgs, function(x){
x$Date<- ymd(ifelse(x$`year(Date)`==2018,
as.character(as.Date(x$`yday(Date)`, origin= "2018-01-01")),
as.character(as.Date(x$`yday(Date)`, origin= "2019-01-01"))))
x$`year(Date)`<- NULL
x$`yday(Date)`<- NULL
return(x)
})
Cor_rain_place<- data.frame(matrix(ncol = 11))
for (i in 1:length(Daiyli_avgs2)) {
Corr1<- cor(Daiyli_avgs2[[i]]$rWRF_acum,
Daiyli_avgs2[[i]]$aport_SMA_diario)
Corr2<- cor(Daiyli_avgs2[[i]]$rWRF_acum,
Daiyli_avgs2[[i]]$aport_mean_diario)
Corr3<- cor(Daiyli_avgs2[[i]]$rWRF_acum,
Daiyli_avgs2[[i]]$aport_mean_SMA_diario)
Corr4<- cor(Daiyli_avgs2[[i]]$rWRF_acum,
Daiyli_avgs2[[i]]$rDHI_acum)
Corr5<- cor(Daiyli_avgs2[[i]]$rWRF_mean,
Daiyli_avgs2[[i]]$aport_SMA_diario)
Corr6<- cor(Daiyli_avgs2[[i]]$rWRF_mean,
Daiyli_avgs2[[i]]$aport_mean_diario)
Corr7<- cor(Daiyli_avgs2[[i]]$rWRF_mean,
Daiyli_avgs2[[i]]$aport_mean_SMA_diario)
Corr8<- cor(Daiyli_avgs2[[i]]$rWRF_mean,
Daiyli_avgs2[[i]]$rDHI_mean)
LON<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][1])
LAT<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][3])
Dist<- distm(c(-7.713948, 42.628577),
c(as.numeric(LON),
as.numeric(LAT)),
fun = distHaversine)
Cor_rain_place[i,]<- as.data.frame(cbind(LON, LAT, round(Corr1, digits = 2),
round(Corr2, digits = 2),
round(Corr3, digits = 2),
round(Corr4, digits = 2),
round(Corr5, digits = 2),
round(Corr6, digits = 2),
round(Corr7, digits = 2),
round(Corr, digits = 2),
Dist))
}
colnames(Cor_rain_place)<- c("LON","LAT","AWRF_aport1",
"AWRF_aport2", "AWRF_aport3",
"AWRF_ADHI","MWRF_aport1",
"MWRF_aport2", "MWRF_aport3",
"MWRF_MDHI","Dist")
Cor_order<- Cor_rain_place[order(Cor_rain_place$Dist),]
Cor_order
Cor_rain_place_ccf<- data.frame(matrix(ncol = 11))
for (i in 1:length(Daiyli_avgs2)) {
ccfr1<- ccf(Daiyli_avgs2[[i]]$rWRF_acum,
Daiyli_avgs2[[i]]$aport_SMA_diario, plot = FALSE)
ccfr2<- ccf(Daiyli_avgs2[[i]]$rWRF_acum,
Daiyli_avgs2[[i]]$aport_mean_diario, plot = FALSE)
ccfr3<- ccf(Daiyli_avgs2[[i]]$rWRF_acum,
Daiyli_avgs2[[i]]$aport_mean_SMA_diario, plot = FALSE)
ccfr4<- ccf(Daiyli_avgs2[[i]]$rWRF_acum,
Daiyli_avgs2[[i]]$rDHI_acum, plot = FALSE)
ccfr5<- ccf(Daiyli_avgs2[[i]]$rWRF_mean,
Daiyli_avgs2[[i]]$aport_SMA_diario, plot = FALSE)
ccfr6<- ccf(Daiyli_avgs2[[i]]$rWRF_mean,
Daiyli_avgs2[[i]]$aport_mean_diario, plot = FALSE)
ccfr7<- ccf(Daiyli_avgs2[[i]]$rWRF_mean,
Daiyli_avgs2[[i]]$aport_mean_SMA_diario, plot = FALSE)
ccfr8<- ccf(Daiyli_avgs2[[i]]$rWRF_mean,
Daiyli_avgs2[[i]]$rDHI_mean, plot = FALSE)
LON<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][1])
LAT<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][3])
Dist<- distm(c(-7.713948, 42.628577),
c(as.numeric(LON),
as.numeric(LAT)),
fun = distHaversine)
Cor_rain_place_ccf[i,]<- as.data.frame(cbind(LON, LAT,
round(max(ccfr1$acf), digits = 2),
round(max(ccfr2$acf), digits = 2),
round(max(ccfr3$acf), digits = 2),
round(max(ccfr4$acf), digits = 2),
round(max(ccfr5$acf), digits = 2),
round(max(ccfr6$acf), digits = 2),
round(max(ccfr7$acf), digits = 2),
round(max(ccfr8$acf), digits = 2), Dist))
print(paste("Mejor corr para un retardo de: ",
ccfr1$lag[which.max(ccfr1$acf)],
ccfr2$lag[which.max(ccfr2$acf)],
ccfr3$lag[which.max(ccfr3$acf)],
ccfr4$lag[which.max(ccfr4$acf)],
ccfr5$lag[which.max(ccfr5$acf)],
ccfr6$lag[which.max(ccfr6$acf)],
ccfr7$lag[which.max(ccfr7$acf)],
ccfr8$lag[which.max(ccfr8$acf)], sep = " "))
}
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -6 -5 -6 0 -6 -5 -6 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -3 -4 -3 0 -3 -4 -3 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -2 0 -3 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -2 0 -3 -2 -2 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de: -13 -12 -13 0 -13 -12 -13 0"
colnames(Cor_rain_place_ccf)<- c("LON","LAT","AWRF_aport1",
"AWRF_aport2", "AWRF_aport3",
"AWRF_ADHI","MWRF_aport1",
"MWRF_aport2","MWRF_aport3",
"MWRF_MDHI","Dist")
Cor_order<- Cor_rain_place_ccf[order(Cor_rain_place_ccf$Dist),]
Cor_order
A continuación extraemos la media y el máximo de cada columna.
#Media de cada columna
colMeans(Cor_rain_place_ccf[3:(length(Cor_rain_place_ccf)-1)])
## AWRF_aport1 AWRF_aport2 AWRF_aport3 AWRF_ADHI MWRF_aport1 MWRF_aport2
## 0.2883721 0.2939535 0.2890698 0.6567442 0.2883721 0.2939535
## MWRF_aport3 MWRF_MDHI
## 0.2890698 0.6567442
#función ColMAX
colMax <- function(data) sapply(data, max, na.rm = TRUE)
colMax(Cor_rain_place_ccf)[3:(length(Cor_rain_place_ccf)-1)]
## AWRF_aport1 AWRF_aport2 AWRF_aport3 AWRF_ADHI MWRF_aport1 MWRF_aport2
## 0.49 0.50 0.49 0.73 0.49 0.50
## MWRF_aport3 MWRF_MDHI
## 0.49 0.73
De aquí podemos extraer las siguientes conclusiones:
#función ColMAX
which_colMax <- function(data) sapply(data, which.max)
which_colMax(Cor_rain_place_ccf)[3:(length(Cor_rain_place_ccf)-1)]
## AWRF_aport1 AWRF_aport2 AWRF_aport3 AWRF_ADHI MWRF_aport1 MWRF_aport2
## 18 18 18 1 18 18
## MWRF_aport3 MWRF_MDHI
## 18 1
Vemos que en cuanto a correlación entre WRF y aportación el mejor punto es el numero 40.
Cor_rain_place_ccf[40,]
Apreciamos que se encuentra a una distancia de 82 km. Esto es peligroso.
En cuanto a la lluvia…
Cor_rain_place_ccf[4,]
Vemos que se encuentra a 52 km.
Cor_rain_place_ccf[order(Cor_rain_place_ccf$AWRF_aport3, decreasing = T),c("AWRF_aport3","AWRF_ADHI", "Dist")]
Por lo visto en la tabla anterior, habrá que buscar un punto que tenga buenas correlaciones para ambas comparaciones y además esté relativamente cerca. En mi opinión estas son las mejores opciones. No sé cual elegir.
Cor_rain_place_ccf[c(19,23,24,14),c("LON","LAT","AWRF_aport3","AWRF_ADHI", "Dist")]
Construimos una lista solo con los puntos las variables que nos interesan. Hemos ido reduciendo el rango de tiempo viendo como las correlaciones aumentaban, ahora que hemos escogido el mejor punto vamos a cojer los datos desde finales de octubre.
clean_data_oct<-lapply(clean_data, function(x){
y<- lapply(x, function(r) {
r<- r[which(r$Date > ymd("2018/10/01")),]
r$prep_hourly<- ifelse(r$prep_hourly < 0,0,r$prep_hourly)
return(r)
})
})
Daiyli_avgs<- lapply(clean_data_oct, function(x){
x[[1]] %>% group_by(year(Date), yday(Date)) %>%
summarize(rWRF_acum= sum(prep_hourly),
rDHI_acum= sum(Lluvia_mm),
rWRF_mean= sum(prep_hourly),
rDHI_mean= sum(Lluvia_mm),
aport_SMA_diario=mean(aport_SMA),
aport_mean_diario=mean(aport_mean),
aport_mean_SMA_diario=mean(aport_mean_SMA),
nivel_SMA_diario=mean(nivel_SMA),
nivel_mean_diario=mean(nivel_mean),
nivel_mean_SMA_diario=mean(nivel_mean_SMA))
})
Daiyli_avgs2<- lapply(Daiyli_avgs, function(x){
x$Date<- ymd(ifelse(x$`year(Date)`==2018,
as.character(as.Date(x$`yday(Date)`, origin= "2018-01-01")),
as.character(as.Date(x$`yday(Date)`, origin= "2019-01-01"))))
x$`year(Date)`<- NULL
x$`yday(Date)`<- NULL
return(x)
})
Daiyli_avgs3<-lapply(Daiyli_avgs2[c(19,23,24,14)], function(x)x[, c("Date","rWRF_acum",
"aport_mean_SMA_diario",
"rDHI_acum")])
head(Daiyli_avgs3[[1]])
range(Daiyli_avgs3[[1]]$Date)
## [1] "2018-10-26" "2019-02-22"
Apreciamos que tenemos un un dato diario entre el 26 de octubre y el 22 de febrero.
lapply(Daiyli_avgs3, function(x, titulo){
k<- (max(c(cumsum(x$rWRF_acum),cumsum(x$rWRF_acum)))/max(c(x$rWRF_acum, x$rDHI_acum)))
ggplot(data=x, aes(x=Date))+
geom_bar(aes(y=rWRF_acum), stat="identity", fill="red", alpha=0.5, size=0.05)+
xlab("Date")+ylab("Lluvia diaria [mm/h]")+theme(panel.background = element_blank(),
panel.grid = element_blank())+
geom_line(aes(y = cumsum(rWRF_acum)/k), group = 1, col="red") +
geom_bar(aes(y=rDHI_acum), stat="identity", fill="blue", alpha=0.4, size=0.05)+
xlab("Date")+ylab("Lluvia diaria [mm/h]")+
theme(panel.background = element_blank(),panel.grid = element_blank())+
geom_line(aes(y = cumsum(rDHI_acum)/k), group = 1, col="blue") +
scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k),
name = " LLuvia acumulada [mm]",
breaks = seq(min(cumsum(x$rDHI_acum)),
max(cumsum(x$rDHI_acum)),
by=50)),
breaks = seq(min(x$rDHI_acum),
max(x$rDHI_acum),
by=10))+
ggtitle(titulo) +
theme(plot.title = element_text(hjust = 0.5))
}, titulo=names(Daiyli_avgs3))
## $`-8.05770874023438__42.6729469299316`
##
## $`-7.32501220703125__42.6962661743164`
##
## $`-7.88290405273438__42.8138198852539`
##
## $`-8.04904174804688__42.5383224487305`
Es interesante apreciar como pese a las altas correlaciones existe una pequeña diferencia al final del periodo en lo que a lluvia acumulada se refiere. Esto nos da pie a investigar cuales son los puntos que calculan mejor la lluvia acumulada.
vec_diff<- as.vector(unlist(lapply(Daiyli_avgs2, function(x){
r<- sum(x$rWRF_acum)-sum(x$rDHI_acum)
return(r)
})))
vec_diff[order(abs(vec_diff))]
## [1] 4.404910 4.729176 -10.388001 -20.714521 31.905542
## [6] -32.582030 36.890532 39.813112 -49.308699 -49.663263
## [11] -55.735768 76.225877 -77.398381 -79.743920 -81.969871
## [16] -87.654145 -91.766186 91.816973 -98.120085 102.573777
## [21] -108.333345 -108.674158 -111.758937 -117.214062 -120.455026
## [26] -127.220242 -130.881619 -135.337830 -136.514876 -147.815890
## [31] 149.886050 -155.839446 -160.733827 -171.574046 180.299144
## [36] -185.876094 191.348111 -193.420657 -213.624562 -220.962254
## [41] -225.314246 -229.023156 -256.119342
Cor_rain_place_ccf[order(abs(vec_diff))[1:10],c("LON","LAT","AWRF_aport3","AWRF_ADHI", "Dist")]
Vemos que entre las dos tablas que hemos sacado… la localización que se repite es la número 14, que se encuentra a 29 km. Esto es muy interesante para realizar luego las predicciones en base a esta localización.
Cor_rain_place_ccf[order(Cor_rain_place_ccf$AWRF_aport3, decreasing = T)[1:10],c("AWRF_aport3","AWRF_ADHI", "Dist")]
Vemos que los puntos que se repiten en ambas listas… buscando mejor correlación con aportación, cercania y mejor resultado con la lluvia acumulada.
best_points<- order(abs(vec_diff))[1:10][order(abs(vec_diff))[1:10]%in%order(Cor_rain_place_ccf$AWRF_aport3, decreasing = T)[1:10]]
best_points
## [1] 28 18 19 32
litros_acum_diff<- vec_diff[best_points]
cbind(litros_acum_diff,
Cor_rain_place_ccf[best_points,c("AWRF_aport3","AWRF_ADHI", "Dist")])
Para un futuro tendremos que elegir entre estos 3 puntos. Añadimos también datos de correlación.
Daiyli_avgs4<-lapply(Daiyli_avgs2[c(14,18,19)], function(x)x[, c("Date","rWRF_acum",
"aport_mean_SMA_diario",
"rDHI_acum")])
lapply(Daiyli_avgs4, function(x, titulo){
k<- (max(c(cumsum(x$rWRF_acum),cumsum(x$rWRF_acum)))/max(c(x$rWRF_acum, x$rDHI_acum)))
ggplot(data=x, aes(x=ymd(Date)))+
geom_bar(aes(y=rWRF_acum), stat="identity", fill="red", alpha=0.5, size=0.05)+
xlab("Date")+ylab("Lluvia diaria [mm/h]")+theme(panel.background = element_blank(),
panel.grid = element_blank())+
geom_line(aes(y = cumsum(rWRF_acum)/k), group = 1, col="red") +
geom_bar(aes(y=rDHI_acum), stat="identity", fill="blue", alpha=0.4, size=0.05)+
xlab("Date")+ylab("Lluvia diaria [mm/h]")+
theme(panel.background = element_blank(),panel.grid = element_blank())+
geom_line(aes(y = cumsum(rDHI_acum)/k), group = 1, col="blue") +
scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k),
name = " LLuvia acumulada [mm]",
breaks = seq(min(cumsum(x$rDHI_acum)),
max(cumsum(x$rDHI_acum)),
by=50)),
breaks = seq(min(x$rDHI_acum),
max(x$rDHI_acum),
by=10))+
ggtitle(titulo) +
theme(plot.title = element_text(hjust = 0.5))+
geom_line(aes(y=aport_mean_SMA_diario/k), col="green",alpha= 0.9)
}, titulo=names(Daiyli_avgs4))
## $`-8.04904174804688__42.5383224487305`
##
## $`-7.31787109375__42.5615844726562`
##
## $`-8.05770874023438__42.6729469299316`
Aquí viene lo gordo. Ya sabemos cuales son los puntos de WRF que mejor correlan con las observaciones aportadas por DHI. Lo primero, será cojer los datos horarios de los puntos anteriormente vistos. Cositas interesantes que quiero probar:
Se me acaba de ocurrir que igual se podría correlar guay la media movil de la precipitación con la aportación. Debido a esa inercia de la que hablamos
Tambien será necesario analizar como se comportar la variación del nivel con la aportación. en principio deberían correlar bien.
path_hist_WRF<- here::here('Data/Parques/Belesar/Historico/Hist_D1D2_DHI_MERGED.RDS')
clean_data<- readRDS(path_hist_WRF)
clean_data_oct<-lapply(clean_data, function(x){
y<- lapply(x, function(r) r[which(r$Date > ymd("2018/10/01")),])
})
fecha_cut<- ymd("2019/01/15")
#cortar en entrenamiento y predicción
cut_train<- lapply(clean_data_oct[c(14,18,19)], function(x,fecha_end){
y<- lapply(x, function(r){
fecha_ini<- ymd("2018/12/01")
Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
return(Jan_data)
})
return(y)
},
fecha_end=fecha_cut)
cut_train<- lapply(cut_train, function(x){
x[[1]][,c("Date","LON", "LAT", "prep_hourly", "Lluvia_mm",
"aport_mean_SMA", "nivel_mean_SMA")]
})
cut_predict<- lapply(clean_data_oct[c(14,18,19)], function(x, fecha_ini){
y<- lapply(x, function(r){
fecha_end<- ymd("2019/02/20")
Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
return(Jan_data)
})
return(y)
},
fecha_ini=fecha_cut)
cut_predict<- lapply(cut_predict, function(x){
x[[1]][,c("Date","LON", "LAT", "prep_hourly", "Lluvia_mm",
"aport_mean_SMA", "nivel_mean_SMA")]
})
Realizamos una gráfica de la diferencia de nivel respecto a la aportación.
x<- cut_train[[1]]
plot(diff(x$nivel_mean_SMA)*20,
ylim = c(-1.5,2),
type = "l",
xlab = "",
xaxt="n",
yaxt="n",
ylab = " ")
lines(x$aport_mean_SMA/100-1, col="red")
Vemos que existe un disparo en la gráfica… eso es porque existe un salto temporal de 3 días en diciembre… en los que el nivel vario 80 cm… cortamos los datos para eliminar ese disparo.
xx<- x[which(x$Date> ymd("2018/12/10")),]
plot(diff(xx$nivel_mean_SMA)*20,
ylim = c(-1.5,2),
type = "l",
xlab = "",
xaxt="n",
yaxt="n",
ylab = " ")
lines(xx$aport_mean_SMA/100-1, col="red")
Comprobamos las cross-correlation de estas dos variables.
ccf_nivel_aport<- ccf(diff(xx$nivel_mean_SMA),
xx$aport_mean_SMA,
lag.max = 1000)
print(paste0("Correlación máxima: ",
max(ccf_nivel_aport$acf), " Se produce para un retraso de ",
ccf_nivel_aport$lag[
which.max(ccf_nivel_aport$acf)], " horas"))
## [1] "Correlación máxima: 0.79527611153213 Se produce para un retraso de -15 horas"
La aportación se ve reflejada en el nivel con un retraso de 15 horas con una correlación bastante alta de casi el 80%.
vec_cor<- vector()
for (i in seq(6, 144, by=2)) {
x<- cut_train[[1]]
x$prep_hourlySMA<- SMA(x$prep_hourly, i)
x<- x[complete.cases(x),]
rain_aport<- ccf(x$prep_hourlySMA,
x$aport_mean_SMA,
lag.max = 1000,
plot = F)
vec_cor[i]<- max(rain_aport$acf)
}
plot(vec_cor,
xlab = "SMA points",
ylab = "Corr")
Se ha decidido que el SMA ha emplear será de 48 puntos. A continuación realizamos una predicción de la aportación empleando diferentes métodos y variables.
path_hist_WRF<- here::here('Data/Parques/Belesar/Historico/Hist_D1D2_DHI_MERGED.RDS')
clean_data<- readRDS(path_hist_WRF)
clean_data_oct<-lapply(clean_data, function(x){
y<- lapply(x, function(r) r[which(r$Date > ymd("2018/10/01")),])
})
fecha_cut<- ymd("2019/01/15")
#cortar en entrenamiento y predicción
cut_train<- lapply(clean_data_oct[c(14,18,19)], function(x,fecha_end){
y<- lapply(x, function(r){
fecha_ini<- ymd("2018/12/01")
Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
return(Jan_data)
})
return(y)
},
fecha_end=fecha_cut)
cut_train<- lapply(cut_train, function(x){
x[[1]][,c("Date","LON", "LAT", "prep_hourly", "Lluvia_mm",
"aport_mean_SMA", "nivel_mean_SMA")]
})
cut_predict<- lapply(clean_data_oct[c(14,18,19)], function(x, fecha_ini){
y<- lapply(x, function(r){
fecha_end<- ymd("2019/02/20")
Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
return(Jan_data)
})
return(y)
},
fecha_ini=fecha_cut)
cut_predict<- lapply(cut_predict, function(x){
x[[1]][,c("Date","LON", "LAT", "prep_hourly", "Lluvia_mm",
"aport_mean_SMA", "nivel_mean_SMA")]
})
x<- cut_train[[1]]
x$aport1<- lag(x$aport_mean_SMA, 24)
x$aport2<- lag(x$aport_mean_SMA, 48)
x$aport3<- lag(x$aport_mean_SMA, 72)
x$aport4<- lag(x$aport_mean_SMA, 96)
x$aport5<- lag(x$aport_mean_SMA, 120)
x$prep_hourly1<- lead(x$prep_hourly, 24)
x$prep_hourlySMA12<- SMA(x$prep_hourly, 12)
x$prep_hourlySMA24<- SMA(x$prep_hourly, 24)
x$prep_hourlySMA36<- SMA(x$prep_hourly, 36)
x$prep_hourlySMA48<- SMA(x$prep_hourly, 48)
x$prep_hourlySMA48_lag<- lag(x$prep_hourlySMA48, 24)
x<- x[complete.cases(x),]
y<- cut_predict[[1]]
y$aport1<- lag(y$aport_mean_SMA, 24)
y$aport2<- lag(y$aport_mean_SMA, 48)
y$aport3<- lag(y$aport_mean_SMA, 72)
y$prep_hourly1<- lead(y$prep_hourly, 24)
y$prep_hourlySMA12<- SMA(y$prep_hourly, 12)
y$prep_hourlySMA24<- SMA(y$prep_hourly, 24)
y$prep_hourlySMA36<- SMA(y$prep_hourly, 36)
y$prep_hourlySMA48<- SMA(y$prep_hourly, 48)
y$prep_hourlySMA48_lag<- lag(y$prep_hourlySMA48, 24)
y<- y[complete.cases(y),]
indexxx<- 285
fecha_ini<- y$Date[indexxx]
fecha_end<- y$Date[indexxx]+ as.difftime(3, units = "days")
y<- y[which(y$Date> fecha_ini & y$Date< fecha_end), ]
fit_1 <- svm(aport_mean_SMA ~ prep_hourlySMA48_lag * prep_hourlySMA48 , data = x)
fit_2 <- svm(aport_mean_SMA ~ prep_hourlySMA48_lag + aport3 + prep_hourlySMA48, data = x)
fit_3 <- svm(aport_mean_SMA ~ prep_hourlySMA48_lag * aport3 * prep_hourlySMA48, data = x)
fit_4 <- lm(aport_mean_SMA ~ prep_hourlySMA48_lag * prep_hourlySMA48, data = x)
fit_5 <- lm(aport_mean_SMA ~ prep_hourlySMA48_lag + aport3 + prep_hourlySMA48, data = x)
fit_6 <- lm(aport_mean_SMA ~ prep_hourlySMA48_lag * aport3 * prep_hourlySMA48, data = x)
uncorrected<-y$aport_mean_SMA
pred_aport1<- predict(fit_1, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
prep_hourlySMA48 =y$prep_hourlySMA48))
pred_aport2<- predict(fit_2, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
aport3= y$aport3,
prep_hourlySMA48 =y$prep_hourlySMA48))
pred_aport3<- predict(fit_3, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
aport3= y$aport3,
prep_hourlySMA48 =y$prep_hourlySMA48))
pred_aport4<- predict(fit_4, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
prep_hourlySMA48 =y$prep_hourlySMA48))
pred_aport5<- predict(fit_5, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
aport3= y$aport3,
prep_hourlySMA48 =y$prep_hourlySMA48))
pred_aport6<- predict(fit_6, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
aport3= y$aport3,
prep_hourlySMA48 =y$prep_hourlySMA48))
plot(uncorrected, x=y$Date, type = "l",
ylim = c(0,600),
xlab = paste0(range(y$Date)),
ylab = "Aportación [m³/s]",
main = "Predicción de aportación \n empleando diferentes métodos")
lines(pred_aport1,x=y$Date, col="red")
lines(pred_aport2,x=y$Date, col="green")
lines(pred_aport3,x=y$Date,col="blue")
lines(pred_aport4,x=y$Date, col="red", lty=2)
lines(pred_aport5,x=y$Date, col="green", lty=2)
lines(pred_aport6,x=y$Date,col="blue", lty=2)
Del anterior gráfico podemos concluir que no vale la pena realizar la predicción aplicando svm (Support Vector Machine), o por lo menos no de esta forma…. Es interesante remarcar que solamente estamos en disposición de realizar un pronóstico a 3 días… dado que el modelo necesita el dato de aportación de hace 3 días para funcionar.
ccf(diff(x$nivel_mean_SMA), x$prep_hourlySMA48, lag.max = 1000)
print(paste0("Correlación máxima: ",
max(ccf_nivel_aport$acf), " Se produce para un retraso de ",
ccf_nivel_aport$lag[
which.max(ccf_nivel_aport$acf)], " horas"))
## [1] "Correlación máxima: 0.79527611153213 Se produce para un retraso de -15 horas"
No tenemos valores de aportación en la página web, solamente nivel… por lo tanto hay que conseguir una buena correlación entre nivel y aportación. O directamente predecir el nivel con la lluvia.
path_hist_WRF<- here::here('Data/Parques/Belesar/Historico/Hist_D1D2_DHI_MERGED.RDS')
clean_data<- readRDS(path_hist_WRF)
clean_data_oct<-lapply(clean_data, function(x){
y<- lapply(x, function(r){
r<- r[which(r$Date > ymd("2018/10/01")),]
r$nivel_mean_SMA1<- c(diff(r$nivel_mean_SMA),0)
return(r)
})
})
fecha_cut<- ymd("2019/01/15")
#cortar en entrenamiento y predicción
cut_train<- lapply(clean_data_oct[c(14,18,19)], function(x,fecha_end){
y<- lapply(x, function(r){
fecha_ini<- ymd("2018/12/01")
Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
return(Jan_data)
})
return(y)
},
fecha_end=fecha_cut)
cut_train<- lapply(cut_train, function(x){
x[[1]][,c("Date","LON", "LAT", "prep_hourly", "Lluvia_mm",
"aport_mean_SMA", "nivel_mean_SMA","nivel_mean_SMA1")]
})
cut_predict<- lapply(clean_data_oct[c(14,18,19)], function(x, fecha_ini){
y<- lapply(x, function(r){
fecha_end<- ymd("2019/02/20")
Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
return(Jan_data)
})
return(y)
},
fecha_ini=fecha_cut)
cut_predict<- lapply(cut_predict, function(x){
x[[1]][,c("Date","LON", "LAT", "prep_hourly", "Lluvia_mm",
"aport_mean_SMA", "nivel_mean_SMA","nivel_mean_SMA1")]
})
x<- cut_train[[1]]
x$nivel1<- lag(x$nivel_mean_SMA1, 24)
x$nivel2<- lag(x$nivel_mean_SMA1, 48)
x$nivel3<- lag(x$nivel_mean_SMA1, 72)
x$nivel4<- lag(x$nivel_mean_SMA1, 96)
x$nivel5<- lag(x$nivel_mean_SMA1, 120)
x$prep_hourly1<- lead(x$prep_hourly, 24)
x$prep_hourlySMA12<- SMA(x$prep_hourly, 12)
x$prep_hourlySMA24<- SMA(x$prep_hourly, 24)
x$prep_hourlySMA36<- SMA(x$prep_hourly, 36)
x$prep_hourlySMA48<- SMA(x$prep_hourly, 48)
x$prep_hourlySMA48_lag<- lag(x$prep_hourlySMA48, 24)
x<- x[complete.cases(x),]
y<- cut_predict[[1]]
y$nivel1<- lag(y$nivel_mean_SMA1, 24)
y$nivel2<- lag(y$nivel_mean_SMA1, 48)
y$nivel3<- lag(y$nivel_mean_SMA1, 72)
y$prep_hourly1<- lead(y$prep_hourly, 24)
y$prep_hourlySMA12<- SMA(y$prep_hourly, 12)
y$prep_hourlySMA24<- SMA(y$prep_hourly, 24)
y$prep_hourlySMA36<- SMA(y$prep_hourly, 36)
y$prep_hourlySMA48<- SMA(y$prep_hourly, 48)
y$prep_hourlySMA48_lag<- lag(y$prep_hourlySMA48, 24)
y<- y[complete.cases(y),]
indexxx<- 285
fecha_ini<- y$Date[indexxx]
fecha_end<- y$Date[indexxx]+ as.difftime(3, units = "days")
y<- y[which(y$Date> fecha_ini & y$Date< fecha_end), ]
fit_1 <- svm(nivel_mean_SMA1 ~ prep_hourlySMA48_lag * prep_hourlySMA48 , data = x)
fit_2 <- svm(nivel_mean_SMA1 ~ prep_hourlySMA48_lag + nivel3 + prep_hourlySMA48, data = x)
fit_3 <- svm(nivel_mean_SMA1 ~ prep_hourlySMA48_lag * nivel3 * prep_hourlySMA48, data = x)
fit_4 <- lm(nivel_mean_SMA1 ~ prep_hourlySMA48_lag * prep_hourlySMA48, data = x)
fit_5 <- lm(nivel_mean_SMA1 ~ prep_hourlySMA48_lag + nivel3 + prep_hourlySMA48, data = x)
fit_6 <- lm(nivel_mean_SMA1 ~ prep_hourlySMA48_lag * nivel3 * prep_hourlySMA48, data = x)
uncorrected<-y$nivel_mean_SMA1
pred_nivel1<- predict(fit_1, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
prep_hourlySMA48 =y$prep_hourlySMA48))
pred_nivel2<- predict(fit_2, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
nivel3= y$nivel3,
prep_hourlySMA48 =y$prep_hourlySMA48))
pred_nivel3<- predict(fit_3, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
nivel3= y$nivel3,
prep_hourlySMA48 =y$prep_hourlySMA48))
pred_nivel4<- predict(fit_4, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
prep_hourlySMA48 =y$prep_hourlySMA48))
pred_nivel5<- predict(fit_5, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
nivel3= y$nivel3,
prep_hourlySMA48 =y$prep_hourlySMA48))
pred_nivel6<- predict(fit_6, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
nivel3= y$nivel3,
prep_hourlySMA48 =y$prep_hourlySMA48))
plot(uncorrected, x=y$Date, type = "l",ylim=c(0, 0.2),
xlab = paste0(range(y$Date)),
ylab = "nivelación [m³/s]",
main = "Predicción de nivelación \n empleando diferentes métodos")
lines(pred_nivel1,x=y$Date, col="red")
lines(pred_nivel2,x=y$Date, col="green")
lines(pred_nivel3,x=y$Date,col="blue")
lines(pred_nivel4,x=y$Date, col="red", lty=2)
lines(pred_nivel5,x=y$Date, col="green", lty=2)
lines(pred_nivel6,x=y$Date,col="blue", lty=2)